home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / splufctr < prev    next >
Text File  |  1994-05-09  |  11KB  |  411 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     Sparse LU factorisation
  29.     See also: sparse.[ch] etc for details about sparse matrices
  30. */
  31.  
  32. #include    <stdio.h>
  33. #include    <math.h>
  34. #include        "sparse2.h"
  35.  
  36.  
  37.  
  38. /* Macro for speedup */
  39. /* #define    sprow_idx2(r,c,hint)    \
  40.    ( ( (hint) >= 0 && (r)->elt[hint].col == (c)) ? hint : sprow_idx((r),(c)) ) */
  41.  
  42.  
  43. /* spLUfactor -- sparse LU factorisation with pivoting
  44.     -- uses partial pivoting and Markowitz criterion
  45.             |a[p][k]| >= alpha * max_i |a[i][k]|
  46.     -- creates fill-in as needed
  47.     -- in situ factorisation */
  48. SPMAT    *spLUfactor(A,px,alpha)
  49. SPMAT    *A;
  50. PERM    *px;
  51. double    alpha;
  52. {
  53.     int    i, best_i, k, idx, len, best_len, m, n;
  54.     SPROW    *r, *r_piv, tmp_row;
  55.     static    SPROW    *merge = (SPROW *)NULL;
  56.     Real    max_val, tmp;
  57.     static VEC    *col_vals=VNULL;
  58.  
  59.     if ( ! A || ! px )
  60.         error(E_NULL,"spLUfctr");
  61.     if ( alpha <= 0.0 || alpha > 1.0 )
  62.         error(E_RANGE,"alpha in spLUfctr");
  63.     if ( px->size <= A->m )
  64.         px = px_resize(px,A->m);
  65.     px_ident(px);
  66.     col_vals = v_resize(col_vals,A->m);
  67.     MEM_STAT_REG(col_vals,TYPE_VEC);
  68.  
  69.     m = A->m;    n = A->n;
  70.     if ( ! A->flag_col )
  71.         sp_col_access(A);
  72.     if ( ! A->flag_diag )
  73.         sp_diag_access(A);
  74.     A->flag_col = A->flag_diag = FALSE;
  75.     if ( ! merge ) {
  76.        merge = sprow_get(20);
  77.        MEM_STAT_REG(merge,TYPE_SPROW);
  78.     }
  79.  
  80.     for ( k = 0; k < n; k++ )
  81.     {
  82.         /* find pivot row/element for partial pivoting */
  83.  
  84.         /* get first row with a non-zero entry in the k-th column */
  85.         max_val = 0.0;
  86.         for ( i = k; i < m; i++ )
  87.         {
  88.         r = &(A->row[i]);
  89.         idx = sprow_idx(r,k);
  90.         if ( idx < 0 )
  91.             tmp = 0.0;
  92.         else
  93.             tmp = r->elt[idx].val;
  94.         if ( fabs(tmp) > max_val )
  95.             max_val = fabs(tmp);
  96.         col_vals->ve[i] = tmp;
  97.         }
  98.  
  99.         if ( max_val == 0.0 )
  100.         continue;
  101.  
  102.         best_len = n+1;    /* only if no possibilities */
  103.         best_i = -1;
  104.         for ( i = k; i < m; i++ )
  105.         {
  106.         tmp = fabs(col_vals->ve[i]);
  107.         if ( tmp == 0.0 )
  108.             continue;
  109.         if ( tmp >= alpha*max_val )
  110.         {
  111.             r = &(A->row[i]);
  112.             idx = sprow_idx(r,k);
  113.             len = (r->len) - idx;
  114.             if ( len < best_len )
  115.             {
  116.             best_len = len;
  117.             best_i = i;
  118.             }
  119.         }
  120.         }
  121.  
  122.         /* swap row #best_i with row #k */
  123.         MEM_COPY(&(A->row[best_i]),&tmp_row,sizeof(SPROW));
  124.         MEM_COPY(&(A->row[k]),&(A->row[best_i]),sizeof(SPROW));
  125.         MEM_COPY(&tmp_row,&(A->row[k]),sizeof(SPROW));
  126.         /* swap col_vals entries */
  127.         tmp = col_vals->ve[best_i];
  128.         col_vals->ve[best_i] = col_vals->ve[k];
  129.         col_vals->ve[k] = tmp;
  130.         px_transp(px,k,best_i);
  131.  
  132.         r_piv = &(A->row[k]);
  133.         for ( i = k+1; i < n; i++ )
  134.         {
  135.         /* compute and set multiplier */
  136.         tmp = col_vals->ve[i]/col_vals->ve[k];
  137.         if ( tmp != 0.0 )
  138.             sp_set_val(A,i,k,tmp);
  139.         else
  140.             continue;
  141.  
  142.         /* perform row operations */
  143.         merge->len = 0;
  144.         r = &(A->row[i]);
  145.         sprow_mltadd(r,r_piv,-tmp,k+1,merge,TYPE_SPROW);
  146.         idx = sprow_idx(r,k+1);
  147.         if ( idx < 0 )
  148.             idx = -(idx+2);
  149.         /* see if r needs expanding */
  150.         if ( r->maxlen < idx + merge->len )
  151.             sprow_xpd(r,idx+merge->len,TYPE_SPMAT);
  152.         r->len = idx+merge->len;
  153.         MEM_COPY((char *)(merge->elt),(char *)&(r->elt[idx]),
  154.             merge->len*sizeof(row_elt));
  155.         }
  156.     }
  157.  
  158.     return A;
  159. }
  160.  
  161. /* spLUsolve -- solve A.x = b using factored matrix A from spLUfactor()
  162.     -- returns x
  163.     -- may not be in-situ */
  164. VEC    *spLUsolve(A,pivot,b,x)
  165. SPMAT    *A;
  166. PERM    *pivot;
  167. VEC    *b, *x;
  168. {
  169.     int    i, idx, len, lim;
  170.     Real    sum, *x_ve;
  171.     SPROW    *r;
  172.     row_elt    *elt;
  173.  
  174.     if ( ! A || ! b )
  175.         error(E_NULL,"spLUsolve");
  176.     if ( (pivot != PNULL && A->m != pivot->size) || A->m != b->dim )
  177.         error(E_SIZES,"spLUsolve");
  178.     if ( ! x || x->dim != A->n )
  179.         x = v_resize(x,A->n);
  180.  
  181.     if ( pivot != PNULL )
  182.         x = px_vec(pivot,b,x);
  183.     else
  184.         x = v_copy(b,x);
  185.  
  186.     x_ve = x->ve;
  187.     lim = min(A->m,A->n);
  188.     for ( i = 0; i < lim; i++ )
  189.     {
  190.         sum = x_ve[i];
  191.         r = &(A->row[i]);
  192.         len = r->len;
  193.         elt = r->elt;
  194.         for ( idx = 0; idx < len && elt->col < i; idx++, elt++ )
  195.         sum -= elt->val*x_ve[elt->col];
  196.         x_ve[i] = sum;
  197.     }
  198.  
  199.     for ( i = lim-1; i >= 0; i-- )
  200.     {
  201.         sum = x_ve[i];
  202.         r = &(A->row[i]);
  203.         len = r->len;
  204.         elt = &(r->elt[len-1]);
  205.         for ( idx = len-1; idx >= 0 && elt->col > i; idx--, elt-- )
  206.         sum -= elt->val*x_ve[elt->col];
  207.         if ( idx < 0 || elt->col != i || elt->val == 0.0 )
  208.         error(E_SING,"spLUsolve");
  209.         x_ve[i] = sum/elt->val;
  210.     }
  211.  
  212.     return x;
  213. }
  214.  
  215. /* spLUTsolve -- solve A.x = b using factored matrix A from spLUfactor()
  216.     -- returns x
  217.     -- may not be in-situ */
  218. VEC    *spLUTsolve(A,pivot,b,x)
  219. SPMAT    *A;
  220. PERM    *pivot;
  221. VEC    *b, *x;
  222. {
  223.     int    i, idx, lim, rownum;
  224.     Real    sum, *tmp_ve;
  225.     /* SPROW    *r; */
  226.     row_elt    *elt;
  227.     static VEC    *tmp=VNULL;
  228.  
  229.     if ( ! A || ! b )
  230.         error(E_NULL,"spLUTsolve");
  231.     if ( (pivot != PNULL && A->m != pivot->size) || A->m != b->dim )
  232.         error(E_SIZES,"spLUTsolve");
  233.     tmp = v_copy(b,tmp);
  234.     MEM_STAT_REG(tmp,TYPE_VEC);
  235.  
  236.     if ( ! A->flag_col )
  237.         sp_col_access(A);
  238.     if ( ! A->flag_diag )
  239.         sp_diag_access(A);
  240.  
  241.     lim = min(A->m,A->n);
  242.     tmp_ve = tmp->ve;
  243.     /* solve U^T.tmp = b */
  244.     for ( i = 0; i < lim; i++ )
  245.     {
  246.         sum = tmp_ve[i];
  247.         rownum = A->start_row[i];
  248.         idx    = A->start_idx[i];
  249.         if ( rownum < 0 || idx < 0 )
  250.         error(E_SING,"spLUTsolve");
  251.         while ( rownum < i && rownum >= 0 && idx >= 0 )
  252.         {
  253.         elt = &(A->row[rownum].elt[idx]);
  254.         sum -= elt->val*tmp_ve[rownum];
  255.         rownum = elt->nxt_row;
  256.         idx    = elt->nxt_idx;
  257.         }
  258.         if ( rownum != i )
  259.         error(E_SING,"spLUTsolve");
  260.         elt = &(A->row[rownum].elt[idx]);
  261.         if ( elt->val == 0.0 )
  262.         error(E_SING,"spLUTsolve");
  263.         tmp_ve[i] = sum/elt->val;
  264.     }
  265.  
  266.     /* now solve L^T.tmp = (old) tmp */
  267.     for ( i = lim-1; i >= 0; i-- )
  268.     {
  269.         sum = tmp_ve[i];
  270.         rownum = i;
  271.         idx    = A->row[rownum].diag;
  272.         if ( idx < 0 )
  273.         error(E_NULL,"spLUTsolve");
  274.         elt = &(A->row[rownum].elt[idx]);
  275.         rownum = elt->nxt_row;
  276.         idx    = elt->nxt_idx;
  277.         while ( rownum < lim && rownum >= 0 && idx >= 0 )
  278.         {
  279.         elt = &(A->row[rownum].elt[idx]);
  280.         sum -= elt->val*tmp_ve[rownum];
  281.         rownum = elt->nxt_row;
  282.         idx    = elt->nxt_idx;
  283.         }
  284.         tmp_ve[i] = sum;
  285.     }
  286.  
  287.     if ( pivot != PNULL )
  288.         x = pxinv_vec(pivot,tmp,x);
  289.     else
  290.         x = v_copy(tmp,x);
  291.  
  292.     return x;
  293. }
  294.  
  295. /* spILUfactor -- sparse modified incomplete LU factorisation with
  296.                         no pivoting
  297.     -- all pivot entries are ensured to be >= alpha in magnitude
  298.     -- setting alpha = 0 gives incomplete LU factorisation
  299.     -- no fill-in is generated
  300.     -- in situ factorisation */
  301. SPMAT    *spILUfactor(A,alpha)
  302. SPMAT    *A;
  303. double    alpha;
  304. {
  305.     int        i, k, idx, idx_piv, m, n, old_idx, old_idx_piv;
  306.     SPROW    *r, *r_piv;
  307.     Real    piv_val, tmp;
  308.     
  309.     /* printf("spILUfactor: entered\n"); */
  310.     if ( ! A )
  311.     error(E_NULL,"spILUfactor");
  312.     if ( alpha < 0.0 )
  313.     error(E_RANGE,"[alpha] in spILUfactor");
  314.     
  315.     m = A->m;    n = A->n;
  316.     sp_diag_access(A);
  317.     sp_col_access(A);
  318.     
  319.     for ( k = 0; k < n; k++ )
  320.     {
  321.     /* printf("spILUfactor(l.%d): checkpoint A: k = %d\n",__LINE__,k); */
  322.     /* printf("spILUfactor(l.%d): A =\n", __LINE__); */
  323.     /* sp_output(A); */
  324.     r_piv = &(A->row[k]);
  325.     idx_piv = r_piv->diag;
  326.     if ( idx_piv < 0 )
  327.     {
  328.         sprow_set_val(r_piv,k,alpha);
  329.         idx_piv = sprow_idx(r_piv,k);
  330.     }
  331.     /* printf("spILUfactor: checkpoint B\n"); */
  332.     if ( idx_piv < 0 )
  333.         error(E_BOUNDS,"spILUfactor");
  334.     old_idx_piv = idx_piv;
  335.     piv_val = r_piv->elt[idx_piv].val;
  336.     /* printf("spILUfactor: checkpoint C\n"); */
  337.     if ( fabs(piv_val) < alpha )
  338.         piv_val = ( piv_val < 0.0 ) ? -alpha : alpha;
  339.     if ( piv_val == 0.0 )    /* alpha == 0.0 too! */
  340.         error(E_SING,"spILUfactor");
  341.  
  342.     /* go to next row with a non-zero in this column */
  343.     i = r_piv->elt[idx_piv].nxt_row;
  344.     old_idx = idx = r_piv->elt[idx_piv].nxt_idx;
  345.     while ( i >= k )
  346.     {
  347.         /* printf("spILUfactor: checkpoint D: i = %d\n",i); */
  348.         /* perform row operations */
  349.         r = &(A->row[i]);
  350.         /* idx = sprow_idx(r,k); */
  351.         /* printf("spLUfactor(l.%d) i = %d, idx = %d\n",
  352.            __LINE__, i, idx); */
  353.         if ( idx < 0 )
  354.         {
  355.         idx = r->elt[old_idx].nxt_idx;
  356.         i = r->elt[old_idx].nxt_row;
  357.         continue;
  358.         }
  359.         /* printf("spILUfactor: checkpoint E\n"); */
  360.         /* compute and set multiplier */
  361.         r->elt[idx].val = tmp = r->elt[idx].val/piv_val;
  362.         /* printf("spILUfactor: piv_val = %g, multiplier = %g\n",
  363.            piv_val, tmp); */
  364.         /* printf("spLUfactor(l.%d) multiplier = %g\n", __LINE__, tmp); */
  365.         if ( tmp == 0.0 )
  366.         {
  367.         idx = r->elt[old_idx].nxt_idx;
  368.         i = r->elt[old_idx].nxt_row;
  369.         continue;
  370.         }
  371.         /* idx = sprow_idx(r,k+1); */
  372.         /* if ( idx < 0 )
  373.         idx = -(idx+2); */
  374.         idx_piv++;    idx++;    /* now look beyond the multiplier entry */
  375.         /* printf("spILUfactor: checkpoint F: idx = %d, idx_piv = %d\n",
  376.            idx, idx_piv); */
  377.         while ( idx_piv < r_piv->len && idx < r->len )
  378.         {
  379.         /* printf("spILUfactor: checkpoint G: idx = %d, idx_piv = %d\n",
  380.                idx, idx_piv); */
  381.         if ( r_piv->elt[idx_piv].col < r->elt[idx].col )
  382.             idx_piv++;
  383.         else if ( r_piv->elt[idx_piv].col > r->elt[idx].col )
  384.             idx++;
  385.         else /* column numbers match */
  386.         {
  387.             /* printf("spILUfactor(l.%d) subtract %g times the ",
  388.                __LINE__, tmp); */
  389.             /* printf("(%d,%d) entry to the (%d,%d) entry\n",
  390.                k, r_piv->elt[idx_piv].col,
  391.                i, r->elt[idx].col); */
  392.             r->elt[idx].val -= tmp*r_piv->elt[idx_piv].val;
  393.             idx++;    idx_piv++;
  394.         }
  395.         }
  396.  
  397.         /* bump to next row with a non-zero in column k */
  398.         /* printf("spILUfactor(l.%d) column = %d, row[%d] =\n",
  399.            __LINE__, r->elt[old_idx].col, i); */
  400.         /* sprow_foutput(stdout,r); */
  401.         i = r->elt[old_idx].nxt_row;
  402.         old_idx = idx = r->elt[old_idx].nxt_idx;
  403.         /* printf("spILUfactor(l.%d) i = %d, idx = %d\n", __LINE__, i, idx); */
  404.         /* and restore idx_piv to index of pivot entry */
  405.         idx_piv = old_idx_piv;
  406.     }
  407.     }
  408.     /* printf("spILUfactor: exiting\n"); */
  409.     return A;
  410. }
  411.